!*********************************************************
!*****   PFS2 v3.0 OZone fire model subroutine library ***
!*********************************************************
!
! VTT Building and Transport                          2003
! VTT Technical Research Centre of Finland
!
! Author:  Timo Korhonen
! Date:    3.9.2003 (version v0.0 started January 10th, 2003)
! Version: PFS2 v2.1 
! Expires: 30 June 2004
!
! DISCLAIMER: This source file comes as it is. There is no waranty
!             that it should operate 'correctly'.
!
!
!*********************************************************
! Version v2.1:
! reads in the OZone '.dat' file generated by the GUI and
! fills the Excel input sheets (the user specified RHR time
! series is not read in, it should be given once again in the
! Excel sheet)
!
!*********************************************************
! Version v2.0:
! Added the BBM/BRE Composite Slab Model
! (This does not change this Fortran file, all calculations
! are done in the Excel macros.)
!*********************************************************
! Version 0.3:
! 26th June 2003: Composite slab model needs file '.cei' for ceiling
!                 temperatures. CallMode = CallMode + 10 if slab
!                 calculation is to be done (iSlabCalc = CallMode/10 and
!                 CallMode = mod(CallMode,10).
!
!*****************************
! Version 0.2:
! 11th April 2003: Changed the arrays: Now Output(nmc,iOutNum)
!                  and RandomNumbers(nmc,iRanNum) (the pfs2 files)
!
! 16th January 2003: save_gastemp-file contains now more data than
!   just the upper layer temperature
!
! 16th January 2003: An expiration date (30 June 2003) is added.
!
! 14th January 2003: A bug is fixed 'SaveSteelTemperatures' is now
!    done better. For each iSteelIter the number of lines is counted
!    and the largest value is used when the save data is written to
!    the disk. (SteelTemps are filled with 1200 C, if there is less
!    than the maximum number of points in the data: SteelPFS.exe stops
!    when steel temperature of 1200 C is reached, so there might be
!    different number of points for different beams.
!*********************************************************
!
! File '.cri' has the transition criteria
!     C1: Tu > T_fo
!     C2: Too hot for the fuel ( T > T_ign )
!     C3: Hot zone too close to the floor
!     C4: A_fire > x % of the floor area
!     BUG in ozonepfs.dll: non-localized fire load, 2 zone calculation
!     and T > Tflashover when Af < Af_crit, then no criteria is printed
!     in the '.cri' file. THIS IS CORRECTED.
!
! File '.op'  has the opening times of windows
!  Opening number, type of variation, criterion value, time of happening
!      2    1    250.00    360.00
!      2    1    350.00    491.02
!      1    1    250.00    684.00
!      1    1    350.00    684.00
! File '.flo' has the floor temperature
! File '.cei' has the ceiling temperature
! File '.WUx' hat the upper wall temperature for wall x
! File '.WLx' hat the lower wall temperature for wall x
!  Time Tgas Tinside ... Toutside
!
!*********************************************************
! Purpose: An interface library (.dll) between Excel and
!          OZone fire zone model.  This interface contains
!          subroutines which:
!            1) writes the .dat input file for OZone (THIS WORKS)
!            1) writes the .da2 input file for SteelTemp (THIS WORKS)
!            2) calls the ozone.exe (THIS WORKS)
!            2) calls the steeltemp.exe (THIS WORKS)
!            3) calls the fireres.exe (THIS WORKS)
!            4) filters OZone's output such that it can be
!               transfered back to Excel (PFS worksheets)
!               (THIS WORKS, but one should decide which numbers are
!                needed in the MC simulations, i.e. which numbers should
!                be filtered back to Excel)
!
!                Now the .pri file information is transfered back
!                and also the steel temperatures (as a function of time)
!                and fire resistance times.
!
!
! This version saves the zone temperature and steel temperature
! data, file names (save_pfs_gastemp.txt and save_pfs_beamtemps.txt)
! CallMode 0: gas temperature
!          1: gas temperature, debug mode
!          2: gas and steel temperatures
!          3: gas and steel temperatures, debug mode
!          4: gas temperatures, save, no debug
!          5: gas and steel temperatures, save, no debug
!          6: gas and steel temperatures, save, debug mode
!
!***********************************************************
!  Calling zone model OZone
!
!  This version writes first a DOS batch file and then runs it.
!  Then it reads the output files and processes them and 
!  returns some variables back to the calling program (Visual
!  Basic macros at Excel)
!
!DEC$ ATTRIBUTES STDCALL,DLLEXPORT :: CallOzone
!DEC$ ATTRIBUTES ALIAS:'CallOzone' :: CallOzone
!***********************************************************
Subroutine CallOzone(cProjName, cRunDir, cPfsDir, cProjectBase, NrO, iVars, CallModeIn, iExtPrint, &
                     iSteelIter, Out, OutSteel, fireres, STMax)

  Use dflib
  Use dfport
  Use dfwin
  Implicit None

!----------------------------------------------
! Argument definitions
!----------------------------------------------
  Character(30) cProjName, cProjectBase
  Character(120) cRunDir, cPfsDir
  Integer(2)    NrO, iVars, CallModeIn, iExtPrint, iSteelIter
  Real(8)       Out(NrO,iVars+1)
  Real(8)       OutSteel(NrO,1+iSteelIter)
  Real(8)       fireres(Max(iSteelIter,1))
  Real(8)       STMax(Max(iSteelIter,1))
!DEC$ ATTRIBUTES REFERENCE :: cProjName, cRunDir, cPfsDir, cProjectBase
!DEC$ ATTRIBUTES REFERENCE :: Out, OutSteel, fireres
!
! Note: Out has an extra column for miscellaneous output
!----------------------------------------------
! Local variables
!----------------------------------------------
  Integer i, j, Iter, iTgasLines
  Integer iDataLines, iLinMax, iSlabCalc, CallMode, iUser
  Integer iochannel, dbgchannel, iostat, ioresultL, savechannel
  Integer i1channel, i2channel, i3channel, sav2channel
  Integer ioresultI, ioresultI4
  Real*8 TimeGasMax
  Character*34 DatFile, cTemp, Da2File, Da3File
  Character*50 cTmpSave
  Character*120 cTmpDir
  Character*120 OzoneDir
  Character*34 MasFile, PriFile, PwrFile, TcfFile
  Character*34 SttFile, NatFile, DbgFile, ResFile
  Character*34 OutFile, BatchFile
  Character*80 line
!  Character*10, Dimension(14) :: cHead, cUnits
  Character*10, Allocatable :: cHead(:), cUnits(:)
  Logical Debug, SaveTemps, Slab, OzoneNotTimeMax
  Real*8, Allocatable :: InData(:,:)
  Character*43 TmpFile, Tmp3File
  Character*39 TmpStr
  Character*8 Buffer
  Character($MAXPATH) dir
!  CHARACTER($MAXPATH) dir
  Integer*4 length
!
! CURRENT DATE
  Integer(2) tmpday, tmpmonth, tmpyear
!
! Process calling variables
!
  integer                     res
  character                   szArgs*512              ! new process arguments buffer & temp pointer
  logical(4)                  fSuccess                ! API return code
  type(T_STARTUPINFO)         si                      ! used for CreateProcess
  type(T_PROCESS_INFORMATION) pi                      ! used for CreateProcess
  type(T_INPUT_RECORD)        ir                      ! used for console input
  type(T_KEY_EVENT_RECORD)    ker

  integer                     dwResult        ! API return code
  integer                     dwCreate        ! new process creation flags
  integer                     SIZEOFSTARTUPINFO
  integer                     SIZESECURITYATTRIBUTES
  integer SIZEOFSzArgs
  character*512 timochar
  character*121 timodir

  SIZEOFSTARTUPINFO = 68
  SIZEOFSzArgs = 512
  SIZESECURITYATTRIBUTES = 12
  dwCreate = CREATE_NEW_CONSOLE
  ! dwCreate = CREATE_NEW_PROCESS_GROUP

!----------------------------------------------
!  Initialize
!----------------------------------------------
  iUser = 0
  CallMode = CallModeIn
  Out = 0.d0
  OutSteel = 0.d0
  fireres = 0.d0
  STMax = 0.d0
  iochannel   = 10
  i1channel   = 11
  i2channel   = 12
  i3channel   = 13
  dbgchannel  = 15
  savechannel = 17
  sav2channel = 18

  ! Pfs_v3p0: Now the location of PfsOzone files is given as a parameter
  OzoneDir = ""
  OzoneDir=cPfsDir(1:index(cPfsDir,char(00))-1)
  OzoneDir=Trim(Ozonedir)

  cTmpDir = ""
  cTmpDir=cRunDir(1:index(cRunDir,char(00))-1)
  cTmpDir=Trim(cTmpDir)

  cTmpSave = ""
  cTmpSave=cProjectBase(1:index(cProjectBase,char(00))-1)
  cTmpSave=Trim(cTmpSave)

  cTemp = ""
  cTemp=cProjName(1:index(cProjName,char(00))-1)
  cTemp=Trim(cTemp)

  DatFile    = Trim(cTemp)//".dat"
  Da2File    = Trim(cTemp)//".da2"
  Da3File    = Trim(cTemp)//".da3"
  PriFile    = Trim(cTemp)//".pri"
  OutFile    = Trim(cTemp)//".OUT"
  MasFile    = Trim(cTemp)//".mas"
  PwrFile    = Trim(cTemp)//".pwr"
  TcfFile    = Trim(cTemp)//".tcf"
  SttFile    = Trim(cTemp)//".stt"
  ResFile    = Trim(cTemp)//".res"
  NatFile    = Trim(cTemp)//".nat"
  DbgFile    = Trim(cTmpSave)//".dbg"
  BatchFile  = Trim(cTemp)//".bat"

  iUser     = CallMode/100
  CallMode  = CallMode - iUser*100
  iUser     = Min(1,iUser)
  ! iUser 0: A/V from catalog
  !       1: A/V user specified
  iSlabCalc = CallMode/10
  iSlabCalc = Min(1,iSlabCalc)    ! Should be 0 or 1
! iSlabCalc 0: Normal steel beam resistance calculation
!           1: BRE composite slab method (Bailey)
  CallMode  = Mod(CallMode,10)    ! Should be 0 - 6
  If ( CallMode .lt. 0 ) CallMode = 0
  If ( CallMode .gt. 6 ) CallMode = 6
! CallMode 0: gas temperature
!          1: gas temperature, debug mode
!          2: gas and steel (and slab) temperatures
!          3: gas and steel (and slab) temperatures, debug mode
!          4: gas temperatures, save, no debug
!          5: gas and steel (and slab) temperatures, save, no debug
!          6: as 5 but debug mode
  If ( CallMode .ge. 4 ) Then
     SaveTemps = .true.
     If ( CallMode .eq. 4 ) CallMode = 0
     If ( CallMode .eq. 5 ) CallMode = 2
     If ( CallMode .ge. 6 ) CallMode = 3
  Else
     SaveTemps = .false.
  End If

  If ((CallMode .eq. 1) .or. (CallMode .eq. 3)) Then
     Debug = .true.
  Else
     Debug = .false.
  End If

  If ( iSlabCalc .eq.1 ) Then
     Slab = .true.
  Else
     Slab = .false.
  End If


!----------------------------------------------
!  Change to Ozone drive and directory
!----------------------------------------------
!  Pfs_v3p0: no more environmental variables
!  ioresultI4  = GETENVQQ ("PFS_v2_DIR", OzoneDir) 
!  OzoneDir = Trim(OzoneDir)
  ioresultL   = CHANGEDIRQQ (Trim(cTmpDir)) 

  Open (dbgchannel, file=Trim(DbgFile), action='readwrite', position='append', form='formatted', iostat=iostat)
  Write (dbgchannel,*) 'Name of the case: ', Trim(cTemp)
!
! CURRENT DATE: Check if this version is expired
!
  Call GETDAT(tmpyear, tmpmonth, tmpday)
  If ( (tmpyear.gt.2008) .or. ( (tmpmonth.gt.6) .and. (tmpyear.gt.2007) ) ) Then
     Write (dbgchannel,*)
     Write (dbgchannel,'(a,i2,''/'',i2.2,''/'',i4.4)') &
          'This version of PFS2 (CallOzone_Pfs_v2p1) has expired on ', &
          6, 30, 2008
     Write (dbgchannel,*)
     Stop
  Else
     If (Debug) Write (dbgchannel,'(a,i2,''/'',i2.2,''/'',i4.4)') &
          'This version of PFS2 (CallOzone_Pfs_v2p1) expires on ', &
          6, 30, 2004
!          tmpmonth, tmpday, tmpyear
  End If

  If (Debug) Then
!     ioresultI = DELFILESQQ (DbgFile) 
!     Open (dbgchannel, file=DbgFile, action='write', form='formatted', iostat=iostat)
     Write (dbgchannel,*)
     Write (dbgchannel,*) '*************************************************'
     Write (dbgchannel,*) 'Name of the case: ', Trim(cTemp)
     Write (dbgchannel,*) '*************************************************'
     Write (dbgchannel,*)
     Write (dbgchannel,*) '**********  CallOzone  ************'

!     If (ioresultI4 .eq. 0 ) Then
!        Write (dbgchannel,*) 'PFS_v2_DIR is not defined'
!     Else If (ioresultI4 .le. 34) Then
     If ( Len_trim(OzoneDir) .gt. 119 ) Then
        Write (dbgchannel,*) 'Too long directory name, maximum is 119 characters'
        Write (dbgchannel,*) Trim(OzoneDir)
     Else
        Write (dbgchannel,*) 'OZonePFS dir is ', Trim(OzoneDir)
!     Else
!        Write (dbgchannel,*) 'PFS_v2_DIR is ', Trim(OzoneDir)
!        Write (dbgchannel,*) 'Error: max 34 characters for PFS_v2_DIR'
     End If

     If ( Len_trim(cTmpDir) .gt. 119 ) Then
        Write (dbgchannel,*) 'Too long directory name, maximum is 119 characters'
        Write (dbgchannel,*) Trim(cTmpDir)
     Else
        Write (dbgchannel,*) 'cTmpDir ioresult is ', ioresultL
        Write (dbgchannel,*) 'Running directory is '
        Write (dbgchannel,*) Trim(cTmpDir)
     End If

     Write (dbgchannel,*) 'Da2File is ', Trim(Da2File)
     Write (dbgchannel,*) 'Da3File is ', Trim(Da3File)
     Write (dbgchannel,*) 'PriFile is ', Trim(PriFile)
     Write (dbgchannel,*) 'OutFile is ', Trim(OutFile)
     Write (dbgchannel,*) 'MasFile is ', Trim(MasFile)
     Write (dbgchannel,*) 'PwrFile is ', Trim(PwrFile)
     Write (dbgchannel,*) 'TcfFile is ', Trim(TcfFile)
     Write (dbgchannel,*) 'SttFile is ', Trim(SttFile)
     Write (dbgchannel,*) 'ResFile is ', Trim(ResFile)
     Write (dbgchannel,*) 'NatFile is ', Trim(NatFile)
     Write (dbgchannel,*) 'DbgFile is ', Trim(DbgFile)

! Number of steel iterations should be less than 30000.
! (If you change this, you must change it also in the WriteSteelData
!  subprogram)
     If ( iSteelIter .gt. 30000 ) Then
        Write (dbgchannel,'( A )') "Too many steel iterations"
        Write (dbgchannel,'( A )') "Only the first 30000 are done correctly"
     End If
  End If

!----------------------------------------
!  Delete existing input and result files
!----------------------------------------
  ioresultI = DELFILESQQ (PriFile) 
  ioresultI = DELFILESQQ (OutFile) 
  If (CallMode .gt. 1 ) Then
     ioresultI = DELFILESQQ(Da2File)
     ioresultI = DELFILESQQ(Da3File)
     ioresultI = DELFILESQQ(SttFile) 
     ioresultI = DELFILESQQ(ResFile) 
     ioresultI = DELFILESQQ(NatFile) 
  End If
  If (iExtPrint .eq. 1 ) Then
     ioresultI = DELFILESQQ (MasFile) 
     ioresultI = DELFILESQQ (PwrFile) 
     ioresultI = DELFILESQQ (TcfFile) 
  End If

   !------------------------------------------------------------------
   !  Execute ozone.exe and steel temp and fire res parts
   !------------------------------------------------------------------
   timochar = Trim(OzoneDir)//"\OzonePFS.exe "//Trim(DatFile)//Char(0)
   timodir = Trim(cTmpDir)//Char(0)
   call ZeroMemory (LOC (si), SIZEOFSTARTUPINFO)
   si%wShowWindow = SW_MINIMIZE
   si%dwFlags = STARTF_USESHOWWINDOW
   fSuccess = .FALSE. ! to enter into while loop
   fSuccess = CreateProcess(NULL_character,    &  ! image file name
                timochar,                        &  ! command line (including program name)
                NULL_security_attributes,      &  ! security for process
                NULL_security_attributes,      &  ! security for main thread
                .FALSE.,                       &  ! new process inherits handles?
                dwCreate,                      &  ! creation flags
                NULL,                          &  ! environment
                timodir,                       &  ! new current DirName
                si,                            &  ! STARTUPINFO structure
                pi)                               ! PROCESSINFORMATION structure

   ! Next does perhaps not work on WindowsXP. It seems that
   ! WindowsXP is not waiting...
   dwResult = WaitForSingleObject(pi%hProcess, -1)
   !
   ! Next might work (from fortran help)
   ! dwResult = WaitForSingleObject(pi%hProcess, INFINITE)

   res = CloseHandle( pi.hThread )
   res = CloseHandle( pi.hProcess )
   If ( Debug ) Then
      Write (dbgchannel,*) 'fSuccess is ',fSuccess
      Write (dbgchannel,*) 'dwResult ',dwResult
      Write (dbgchannel,'("ozonepfs.exe executed")') 
   End If

  Do Iter = 1, iSteelIter
     If ( Iter .le. 30000 ) Then
        Write (Buffer, '(i8)') Iter
        TmpFile = Trim(cTemp)//"-"//Trim(Adjustl(Buffer))
        TmpStr  = Trim(TmpFile)
        TmpFile = Trim(TmpStr)//".da2"
        Tmp3File = Trim(TmpStr)//".da3"
     Else
        TmpFile = Trim(cTemp)//"-ERROR"//".da2"
        Tmp3File = Trim(cTemp)//"-ERROR"//".da3"
        TmpStr  = Trim(cTemp)//"-ERROR"
     End If
     ioresultI = DELFILESQQ(Trim(Da3File))
     ioresultI = DELFILESQQ(Trim(Da2File))
     ioresultI = RENAMEFILEQQ(Trim(TmpFile), Trim(Da2File))
     ioresultI = RENAMEFILEQQ(Trim(Tmp3File), Trim(Da3File))
     timochar = ""
     ! Note: SteelPFS.exe does not work with user defined A/V. Use the
     ! older version of the steel routine.
     If ( iUser .eq. 1 ) Then
        timochar = Trim(OzoneDir)//"\SteelPFS_user.exe "//Trim(Da2File)//Char(0)
     Else
        timochar = Trim(OzoneDir)//"\SteelPFS.exe "//Trim(Da2File)//Char(0)
     End If
     If (Debug) Then
        Write (dbgchannel,*) 'timochar ',timochar
        Write (dbgchannel,*) 'cTemp ',cTemp
        Write (dbgchannel,*) 'TmpFile ',Tmpfile
        Write (dbgchannel,*) 'Tmp3File ',Tmp3file
        Write (dbgchannel,*) 'Da2File ',Da2file
        Write (dbgchannel,*) 'Da3File ',Da3file
     End If

     call ZeroMemory (LOC (si), SIZEOFSTARTUPINFO)
     si%dwFlags = STARTF_USESHOWWINDOW
     si%wShowWindow = SW_MINIMIZE
     fSuccess = .FALSE. ! to enter into while loop
     fSuccess = CreateProcess(null_character,    &  ! image file name
          timochar,                        &  ! command line (including program name)
          NULL_security_attributes,      &  ! security for process
          NULL_security_attributes,      &  ! security for main thread
          .FALSE.,                       &  ! new process inherits handles?
          dwCreate,                      &  ! creation flags
          NULL,                          &  ! environment
          timodir,                      &  ! new current DirName
          si,                            &  ! STARTUPINFO structure
          pi)                               ! PROCESSINFORMATION structure

     dwResult = WaitForSingleObject(pi%hProcess, -1)
     res = CloseHandle( pi.hThread )
     res = CloseHandle( pi.hProcess )
     If ( Debug ) Then
        Write (dbgchannel,*) 'fSuccess is ',fSuccess
        Write (dbgchannel,*) 'dwResult ',dwResult
        Write (dbgchannel,'("SteelPFS.exe executed")') 
     End If

     timochar = ""
     timochar = Trim(OzoneDir)//"\FireResPFS.exe "//Trim(Da3File)// Char(0)
     call ZeroMemory (LOC (si), SIZEOFSTARTUPINFO)
     si%dwFlags = STARTF_USESHOWWINDOW
     si%wShowWindow = SW_MINIMIZE
     fSuccess = .FALSE. ! to enter into while loop
     fSuccess = CreateProcess(null_character,    &  ! image file name
          timochar,                        &  ! command line (including program name)
          NULL_security_attributes,      &  ! security for process
          NULL_security_attributes,      &  ! security for main thread
          .FALSE.,                       &  ! new process inherits handles?
          dwCreate,                      &  ! creation flags
          NULL,                          &  ! environment
          timodir,                      &  ! new current DirName
          si,                            &  ! STARTUPINFO structure
          pi)                               ! PROCESSINFORMATION structure
     
     dwResult = WaitForSingleObject(pi%hProcess, -1)
     res = CloseHandle( pi.hThread )
     res = CloseHandle( pi.hProcess )
     If ( Debug ) Then
        Write (dbgchannel,*) 'fSuccess is ',fSuccess
        Write (dbgchannel,*) 'dwResult ',dwResult
        Write (dbgchannel,'("FireResPFS.exe executed")') 
     End If
     ioresultI = DELFILESQQ(Trim(TmpStr)//".stt")
     ioresultI = DELFILESQQ(Trim(TmpStr)//".res")
     ioresultI = RENAMEFILEQQ(Trim(SttFile), Trim(TmpStr)//".stt")
     ioresultI = RENAMEFILEQQ(Trim(ResFile), Trim(TmpStr)//".res")
  End Do

!--------------------------------------------
!  Open and Read the result file
!--------------------------------------------
!  If (Debug) Then
!     ioresultI = DELFILESQQ (DbgFile) 
!     Open (dbgchannel, file=DbgFile, action='readwrite', position='append', form='formatted', iostat=iostat)
!  End If

  Open (iochannel, file=Trim(PriFile), action='read', form='formatted', iostat=iostat)

! Calculate the number of lines in the output
! The 'NrO' parameter is the maximum number of lines which
! is used in the Excel sheet. (It is now 100, which should
! be enough to plot nice figures.)
  i = 0
  Do While (.not. EOF(iochannel))
     i = i + 1
     Read (iochannel, *)
  End Do
  Rewind (iochannel)
  ! First line is variable names, second line units
  iDataLines = i - 2
  If (Debug) Then
     Write (dbgchannel,*) 'Number of Pri Data Lines and NrO are ', iDataLines, NrO
  End If
  Allocate(cHead(iVars+1))
  Allocate(cUnits(iVars+1))
  Allocate(InData(iDataLines,iVars+1))
  InData = 0.d0
  Call ReadOzoneOutput(Debug, iochannel, dbgchannel, &
       InData, cHead, cUnits, iDataLines, iVars)

  OzoneNotTimeMax = .False.
  iTgasLines = NrO
  If ( iDataLines .lt. NrO ) Then 
     iTgasLines = iDataLines
     OzoneNotTimeMax = .True.
  End If
  TimeGasMax = InData(iTgasLines,1)

  Do i = 1, Min0(iDataLines,NrO)
     Out(i,:) = InData(i,:)
  End Do
  Deallocate(InData)
  Close(iochannel)
  If ( Debug ) Then
     Write (dbgchannel,*) 'PriFile closed'
  End If
!
! Write gas temperatures at the save file. Each iteration data is written
! at the end of the existing file (so remember delete the file before the
! first iteration).
!
  If ( SaveTemps ) Then
     Open (savechannel, file=Trim(cTmpSave)//'_gastemp.txt', action='readwrite', position='append', form='formatted', iostat=iostat)
     Write (savechannel,*) 'Name of the iteration: ', Trim(cTemp)
     Write (savechannel,*) iDataLines
     Write (savechannel, fmt='(a8,13a10)') (cHead(i), i = 1,iVars)
     Write (savechannel, fmt='(a8,13a10)') (cUnits(i), i = 1,iVars)
     Do i = 1, Min0(iDataLines,NrO)
        Write (savechannel,fmt='(f8.1,5f10.3,2f10.0,6f10.3)') (Out(i,j), j = 1,iVars)
     End Do
     Close (savechannel)
  End If
  Deallocate(cUnits)
  Deallocate(cHead)
!
! Steel temperatures and fire resistance
!  
  Do Iter = 1, iSteelIter
     If ( Iter .le. 30000 ) Then
        Write (Buffer, '(i8)') Iter
        TmpFile = Trim(cTemp)//"-"//Trim(Adjustl(Buffer))//".stt"
     Else
        TmpFile = Trim(cTemp)//"-ERROR"//".stt"
     End If

     If (Debug) Then
        Write(dbgchannel,*) 'SteelIter ',Iter
        Write(dbgchannel,*) 'TmpFile Stt ',TmpFile
     End If
     Open (i1channel, file=Trim(TmpFile), action='read', form='formatted', iostat=iostat)

     iLinMax = 1
     If ( Iter .eq. 1 ) Then
        OutSteel = 0.d0
     End If
! Calculate the number of lines in the output
! The 'NrO' parameter is the maximum number of lines which
! is used in the Excel sheet. (It is now 100, which should
! be enough to plot nice figures.)
     i = 0
     Do While (.not. EOF(i1channel))
        i = i + 1
        Read (i1channel, *)
     End Do
     Rewind (i1channel)
     ! First line is title, second line names and units
     iDataLines = i - 2
     If ( iDataLines .le. NrO ) Then
        If (Debug) Then
           Write (dbgchannel,*) 'Number of Stt Data Lines is ', iDataLines
           Write (dbgchannel,*) 'Number of Pri Data Lines is ', iTgasLines, TimeGasMax
        End If
     Else
        iDataLines = NrO
        If (Debug) Then
           Write (dbgchannel,*) 'Number of Stt Data Lines is ', iDataLines
           Write (dbgchannel,*) 'Number of Pri Data Lines is ', iTgasLines, TimeGasMax
           Write (dbgchannel,*) 'Number of rows in Excel is ', NrO
           Write (dbgchannel,*) 'iDataLines IS SET EQUAL TO ', NrO
        End If
     End If
     iLinMax = Max(iLinMax,iDataLines)

     Allocate(InData(iDataLines,3))
     Call ReadSteelOutput(Debug, i1channel, dbgchannel, InData, &
                          iDataLines, STMax(Iter))
     Do i = 1, Min0(iDataLines,NrO)
!        If (Iter .eq. 1) Then
           OutSteel(i,1) = InData(i,1)
!        End If
        OutSteel(i,1+Iter) = InData(i,3)
     End Do
     ! The steel calculation stops when a steel temperature of 1200 C
     ! is reached, so fill the array with 1200 C. Fill also time axis.
     Do i = Min0(iDataLines,NrO) + 1, NrO
!        OutSteel(i,1+Iter) = 1200.d0
        OutSteel(i,1+Iter) = OutSteel(iDataLines,1+Iter)
        OutSteel(i,1) = OutSteel(iDataLines,1) + (i-iDataLines)* &
             ( OutSteel(iDataLines,1)-OutSteel(iDataLines-1,1) )
     End Do

     Deallocate(InData)
     Close(i1channel)
     If ( Debug ) Then
        Write (dbgchannel,*) 'After ReadSteelOutput'
        Write (dbgchannel,*) Iter, STMax(Iter)
        Write (dbgchannel,*) 'SttFile closed'
        ioresultI4 = DELFILESQQ(TmpFile)
     Else
        ioresultI4 = DELFILESQQ(TmpFile)
     End If

     ! Open the fire resistance file
     If ( Iter .le. 30000 ) Then
        Write (Buffer, '(i8)') Iter
        TmpFile = Trim(cTemp)//"-"//Trim(Adjustl(Buffer))//".res"
     Else
        TmpFile = Trim(cTemp)//"-ERROR"//".res"
     End If
     If (Debug) Then
        Write(dbgchannel,*) 'SteelIter ',Iter
        Write(dbgchannel,*) 'TmpFile Res ',TmpFile
     End If

     Open (i2channel, file=Trim(TmpFile), action='read', form='formatted', iostat=iostat)
     Read (i2channel, '(A)') line
     Read (i2channel, *) fireres(Iter)

     Close (i2channel)
     If ( Debug ) Then
        Write (dbgchannel,*) 'FireResTime is ',fireres(Iter)
        Write (dbgchannel,*) 'ResFile closed'
        ioresultI4 = DELFILESQQ(TmpFile)
     Else
        ioresultI4 = DELFILESQQ(TmpFile)
     End If

  End Do  ! iSteelIter loop ends
!
! Write steel temperatures at the save file. Each iteration data is written
! at the end of the existing file (so remember delete the file before the
! first iteration).
!
  If ( iSteelIter .ge. 1 ) Then
     If ( SaveTemps ) Then
        Open (sav2channel, file=Trim(cTmpSave)//'_steeltemp.txt', action='readwrite', position='append', form='formatted', iostat=iostat)
        Write (sav2channel,*) 'Name of the iteration: ', Trim(cTemp)
        Write (sav2channel,Fmt='(3i10,f10.1)') iLinMax, iSteelIter, iTgasLines, TimeGasMax
!        Write (sav2channel,*) iDataLines
        Write (sav2channel,*) 'Time(s)  Tbeam1(C) Tbeam2(C)   etc'
!        Do i = 1, Min0(iLinMax,NrO)
        ! Maximun time is at most the maximum time in the gas temperature calculation
        ! (If Ozone crashes during calculation, do not calculate steel temperature after
        ! the crash.)
        Do i = 1, Min0(iLinMax,NrO)
           Write (sav2channel,fmt='(f8.1,10f10.3)') OutSteel(i,1), &
                (OutSteel(i,1+Iter), Iter=1,iSteelIter)
        End Do
        Close (sav2channel)
     End If
  End If
!------------------------------------------------------------------
!  Close debugfile
!------------------------------------------------------------------
  If ( Debug ) Then
     Write (dbgchannel,*) 'DbgFile closed'
     Write (dbgchannel,*) '***** CallOzone ends *****'
  Else
     ioresultI4 = DELFILESQQ(DatFile)
     ioresultI4 = DELFILESQQ(PriFile)
     ioresultI4 = DELFILESQQ(OutFile)
     If (CallMode .gt. 1 ) Then
        ioresultI4 = DELFILESQQ(Da2File)
        ioresultI4 = DELFILESQQ(Da3File)
        ioresultI4 = DELFILESQQ(SttFile) 
        ioresultI4 = DELFILESQQ(ResFile) 
        ioresultI4 = DELFILESQQ(NatFile) 
     End If
     If (iExtPrint .eq. 1 ) Then
        ioresultI = DELFILESQQ (MasFile) 
        ioresultI = DELFILESQQ (PwrFile) 
        ioresultI = DELFILESQQ (TcfFile) 
     End If
  End If
  Close (dbgchannel)
  Return

Contains

!***********************************************************
! Read the output files of OZone and processes them also
!***********************************************************
  Subroutine ReadOzoneOutput( Debug, iochannel, dbgchannel, &
       InData, cHead, cUnits, iDataLines, iVars)

    Implicit None
!----------------------------------------------
! Passed variables
!----------------------------------------------
    Integer iochannel, dbgchannel, iVars, iDataLines
    Logical Debug
    Real*8 InData(iDatalines,iVars+1)
    Character*10 cHead(iVars+1), cUnits(iVars+1)
!----------------------------------------------
! Local variables
!----------------------------------------------
    Integer i, j
    Real*8 Zi1, Zi2, Tup1, Tup2, TempVal, TempTime
!----------------------------------------------
!  Initialize
!----------------------------------------------

    Read (iochannel, fmt='(a8,13a10)') (cHead(i), i = 1,iVars)
    Read (iochannel, fmt='(a8,13a10)') (cUnits(i), i = 1,iVars)
    If (Debug) Then
       Write (dbgchannel,*) '**** ReadOzoneOutput *****'
       Write (dbgchannel,'(a8,13a10)') (cHead(i), i = 1,iVars)
       Write (dbgchannel,'(a8,13a10)') (cUnits(i), i = 1,iVars)
    End If
    Do i = 1, iDataLines
       Read (iochannel,'(f8.1,5f10.3,2f10.0,6f10.3)') (InData(i,j), j=1,iVars)
    End Do
!    If (Debug) Then
!       Write (6,*) 'iVars is ',iVars
!       Do i = 1, 2
!          Write (dbgchannel,'(f8.1,5f10.3,2f10.0,6f10.3)') (InData(i,j), j=1,iVars)
!       End Do
!    End If
    Return
  End Subroutine ReadOzoneOutput

!***********************************************************
! Read the output files of SteelTemp and processes them also
!***********************************************************
  Subroutine ReadSteelOutput(Debug, iochannel, dbgchannel, &
       InData, iDataLines, STMax)

    Implicit None
!----------------------------------------------
! Passed variables
!----------------------------------------------
    Integer iochannel, dbgchannel, iDataLines
    Logical Debug
    Real*8 STMax
    Real*8 InData(iDatalines,3)
!----------------------------------------------
! Local variables
!----------------------------------------------
    Integer i, j
    Character*10 cTmp
    Real*8 Zi1, Zi2, Tup1, Tup2, TempVal, TempTime
!----------------------------------------------
!  Initialize
!----------------------------------------------
    InData = 0.d0
    Read (iochannel, fmt='(a)') cTmp
    Read (iochannel, fmt='(a)') cTmp
    Do i = 1, iDataLines
       Read (iochannel,'(f8.1,2f10.2)') (InData(i,j), j=1,3)
    End Do
!    If (Debug) Then
!       Do i = 1, iDataLines
!          Write (dbgchannel,'(f8.1,2f10.2)') (InData(i,j), j=1,3)
!       End Do
!    End If
    STMax = Maxval(InData(:,3))
    Return
  End Subroutine ReadSteelOutput

End Subroutine CallOzone
!***********************************************************

!***********************************************************
!
!  Write OZone Dat -file
!
!  26th September 2002
!  This works if:
!     -  t2-fire is used
!     -  exactly four walls
!     -  each partition has 7 elements
!     -  max. number of vents is 10
!     -  max. number of materials is 10
!
!DEC$ ATTRIBUTES STDCALL,DLLEXPORT        :: WriteOzoneData
!DEC$ ATTRIBUTES ALIAS : "WriteOzoneData"   :: WriteOzoneData
!***********************************************************
Subroutine WriteOzoneData (cProjName, cRunDir, iStrategy,  &
     TFlashOver, rTFuelIgni, rPercAfionAf, rPercZsonH, &
     rAfloor, rHeight, cHeight, zFireElev, rfuelheight, &
     Tini, rPinitial, rDHc, repsilglass, &
     xtimestep, printstep, timemax, &
     iceiling, ifloor, nbwall, nbmat, nbopen, irmf, &
     nbptrmf, iceiltype, icomb, iairent, nbsmokeext, &
     nbHVgroup, iExtPrint, &
     ropper1, ropT2, ropper2, ropT3, ropper3, &
     roptime2, ropT0, openings, materials, fvents, hvents, &
     MaxWallNodes, nbeleC, nbeleF, walldata, partdata, rmfmax, tplateau, &
     tendplate, tendfire, rfirearea, RHRf, fireseries )

  Use dflib
  Use dfport
  Implicit None
!----------------------------------------------
! Arguments
!----------------------------------------------
  Character(30)  cProjName
  Character(120) cRunDir
  Integer(2)     iStrategy
  Real(8)        TFlashOver, rTFuelIgni, rPercAfionAf, rPercZsonH  
  Real(8)        rAfloor, rHeight, cHeight, zFireElev, rfuelheight
  Real(8)        Tini, rPinitial, rDHc, repsilglass
  Real(8)        xtimestep, printstep, timemax
  Integer(2)     iceiling, ifloor, nbwall, nbmat, nbopen, irmf
  Integer(2)     nbptrmf, iceiltype, icomb, iairent, nbsmokeext
  Integer(2)     nbHVgroup, iExtPrint
  Real(8)        ropper1(3), ropT2(2), ropper2(2), ropT3(2), ropper3(2)
  Real(8)        roptime2, ropT0
  Real(8)        openings(10,7), materials(10,7), fvents(3,5), hvents(3,4)
  Integer(2)     MaxWallNodes, nbeleC, nbeleF
  Real(8)        walldata(2*nbwall)
  Real(8)        partdata(MaxWallNodes,2*nbwall+4)
  Real(8)        rmfmax, tplateau, tendplate, tendfire, rfirearea, RHRf
  Real(8)        fireseries(nbptrmf,4)
!DEC$ ATTRIBUTES REFERENCE :: cProjName, cRunDir
!DEC$ ATTRIBUTES REFERENCE :: ropper1, ropT2, ropper2, ropT3, ropper3
!DEC$ ATTRIBUTES REFERENCE :: openings, materials, fvents, hvents
!DEC$ ATTRIBUTES REFERENCE :: walldata, partdata, fireseries
!----------------------------------------------
! Local variables
!----------------------------------------------
  Integer       iochannel, iostat, ioresultL
  Character*34  DatFile
  Character*30  cTemp
  Character*120 cTmpDir
  Integer       i, j, i1, nbeleW
!
! CURRENT DATE
  Integer(2) tmpday, tmpmonth, tmpyear
!

  iochannel = 12
  cTmpDir = ""
  cTmpDir=cRunDir(1:index(cRunDir,char(00))-1)
  cTmpDir=Trim(cTmpDir)
  cTemp = ""
  cTemp=cProjName(1:index(cProjName,char(00))-1)
  cTemp=Trim(cTemp)
  DatFile   = Trim(cTemp)//".dat"
!
!----------------------------------------------
!  Change to the run directory
!----------------------------------------------
  ioresultL   = CHANGEDIRQQ (Trim(cTmpDir)) 
!
!----------------------------------------------
!  Create datfile
!----------------------------------------------
!
  Open (iochannel, file = Trim(DatFile), action = 'write', form='formatted',iostat = iostat)
!
! CURRENT DATE: Check if this version is expired
!
  Call GETDAT(tmpyear, tmpmonth, tmpday)
  If ( (tmpyear.gt.2008) .or. ( (tmpmonth.gt.6) .and. (tmpyear.gt.2007) ) ) Then
     Write (iochannel,'(a,i2,''/'',i2.2,''/'',i4.4)') &
          'This version of PFS2 (CallOzone_Pfs_v2p1) has expired on ', &
          6, 30, 2008
     Stop
  End If

  Write (iochannel,'(A)') Trim(cTemp)
  Write (iochannel,'(" iStrategy")')
  Write (iochannel,'(i10)') iStrategy
  Write (iochannel,'("     TFlashOver     rTFuelIgni   rPercAfionAf     rPercZsonH ")')
  Write (iochannel,'(4f15.2)') TFlashOver, rTFuelIgni, rPercAfionAf, rPercZsonH
  Write (iochannel,'("   rAfloor        rHeight      zFireElev    rfuelheight ")')
  Write (iochannel,'(f10.2,3f15.2)') rAfloor, rHeight, zFireElev, rfuelheight
  Write (iochannel,'("      Tini      rPinitial           rDHc    repsilglass ")')
  Write (iochannel,'(f10.2,3f15.2)') Tini, rPinitial, rDHc, repsilglass
  Write (iochannel,'(" xtimestep      printstep        timemax ")')
  Write (iochannel,'(f10.2,2f15.2)') xtimestep, printstep, timemax

  ! Switches and numbers of different objects
  Write (iochannel,'(a)') '  iceiling    ifloor    nbwall     nbmat    nbopen' // &
       '      irmf   nbptrmf iceiltype     icomb   iairent ' // &
       'nbsmokeext nbHVgroup iExtPrint '
  If ( irmf .eq. 0 ) Then
     Write (iochannel,'(13i10)') iceiling, ifloor, nbwall, nbmat, nbopen, irmf, 0, &
          iceiltype, icomb, iairent, nbsmokeext, nbHVgroup, iExtPrint
  Else
     Write (iochannel,'(13i10)') iceiling, ifloor, nbwall, nbmat, nbopen, irmf, nbptrmf, &
          iceiltype, icomb, iairent, nbsmokeext, nbHVgroup, iExtPrint
  End If

  ! Finite element data for the partitions
  Write (iochannel,'(a)') '   CEILING'
  Write (iochannel,'(a)') '    nbeleC'
  Write (iochannel,'(a)') '    rlelem   nmat'
  If (iceiltype .eq. 0) Then
     Write (iochannel,'(i10)') nbeleC
  Else
     Write (iochannel,'(i10,f7.1)') nbeleC,cHeight
  End If
  Do j = 1, nbeleC
     Write (iochannel,'(f10.6,i7)') partdata(j,1), Nint(partdata(j,2))
  End Do

  Write (iochannel,'(a)') '     FLOOR'
  Write (iochannel,'(a)') '    nbeleF'
  Write (iochannel,'(a)') '    rlelem   nmat'
  If (iceiltype .eq. 0 .or. iceiltype .eq. 2) Then
     Write (iochannel,'(i10)') nbeleF
  Else
     Write (iochannel,'(i10,f7.1)') nbeleF,cHeight
  End If
  Do j = 1, nbeleF
     Write (iochannel,'(f10.6,i7)') partdata(j,3),  Nint(partdata(j,4))
  End Do

  Write (iochannel,'(a)') '     WALLS'
  Write (iochannel,'(a)') 'walllength nbeleW'
  Write (iochannel,'(a)') '    rlelem   nmat'
  Do i = 1, nbwall
     i1 = 2*(i-1) + 1 
     nbeleW = Nint(walldata(i1+1))
     Write (iochannel,'(f10.2,i7)') walldata(i1), nbeleW
     i1 = i1 + 4 ! now the first 4 points are used for floor and ceiling
     Do j = 1, nbeleW
        Write (iochannel,'(f10.6,i7)') partdata(j,i1),  Nint(partdata(j,i1+1))
     End Do
  End Do

  ! Materials used for the walls, floor and ceiling
  Write (iochannel,'(a)') 'MATERIAL FOR PARTITIONS'
  Write (iochannel,'(a)') '       rho        rk        rc repsilint' // &
       ' repsilext     rhint     rhext'
  Do i = 1, nbmat
     Write (iochannel,'(f10.1,f10.3,f10.1,4f10.3)') (materials(i,j), j=1,7)
  End Do

  ! Openings (doors and windows)
  Write (iochannel,'("     rsill   rsoffit    rwidth    rcvent   nopwall   iopevol     iadia ")')
  Do i = 1, nbopen
     Write (iochannel,'(7f10.2)') (openings(i,j), j=1,7)
  End Do

  ! RHR data
  If ( irmf .eq. 0 ) Then
     Write (iochannel,'("    rmfmax  tplateau tendplate  tendfire rfirearea      RHRf ")')
     Write (iochannel,'(f10.7,3f10.2,f10.3,f10.0)') rmfmax, tplateau, tendplate, &
          tendfire, rfirearea, 1000*RHRf
  Else
     Write (iochannel,'("      time       RHR       rmf    rfArea ")')
     Do i = 1, nbptrmf
        Write (iochannel,'(i10,i15,f15.3,f15.2)') Int(fireseries(i,1)),&
        Int(fireseries(i,2)),fireseries(i,3),fireseries(i,4)
     End Do
  End If

  ! Openings: Opening criteria:
  ! first line:  temperature dependent, stairs
  ! second line: temperature dependent, piecewise linear
  ! third line:  time dependent
  ! fourth line: temperature dependent, single step, i.e. closed vs open
  Write (iochannel,'("   ropper1     ropT2   ropper2     ropT3   ropper3 ")')
  Write (iochannel,'(5f10.2)') ropper1(1), ropT2(1), ropper2(1), ropT3(1), ropper3(1)
  Write (iochannel,'(5f10.2)') ropper1(2), ropT2(2), ropper2(2), ropT3(2), ropper3(2)
  Write (iochannel,'("   ropper1  roptime2 ")')
  Write (iochannel,'(2f10.2)') ropper1(3), roptime2
  Write (iochannel,'("               ropT0 ")')
  Write (iochannel,'("          " f10.2)') ropT0

  ! Forced vents (vertical)
  Write (iochannel,'("       zFV    DiamFV       VFV   signMFV iopevolFV ")')
  Do i = 1, nbsmokeext
     Write (iochannel,'(5f10.2)') (fvents(i,j), j=1,5)
  End Do

  ! Horizontal vents (smoke hatches)
  Write (iochannel,'("      nbHV   rAreaHV iopevolHV ")')
  Do i = 1, nbHVgroup
     Write (iochannel,'(5f10.2)') (hvents(i,j), j=1,4)
  End Do
  Write (iochannel,*)
!--------------------------------------------------
! Ready to finish
!--------------------------------------------------
  Close (iochannel)
!--------------------------------------------------
! Finished so return
!--------------------------------------------------
  Return
End Subroutine WriteOzoneData
!***********************************************************
! Timo Korhonen/VTT
! 6th November, 2002
!
! This routine writes the input data file '.da2' for the
! steel temperature calculation. Works for the unprotected
! and protected cases.
!
!DEC$ ATTRIBUTES STDCALL,DLLEXPORT        :: WriteSteelData
!DEC$ ATTRIBUTES ALIAS : "WriteSteelData"   :: WriteSteelData
!***********************************************************
Subroutine WriteSteelData (cProjName, cRunDir, sec_name, Iter, Heat_mode, &
                           Protection, rhoa, &
                           rHbeam, eps1res, Timemax, Timestep, rpasimp, &
                           Hbeam, Rbeam, AmsurV, Rdp, ProtRows, &
                           Analysis, GammaMFi, Steel_grade, NFid, Lfy, &
                           Lfz, K_index, K1_bend, K2_bend, Lbraced, &
                           LoadType, Psi, Diagram, Zg, User, Rolled, &
                           B_user, Tf_user, Hi_user, Tw_user, &
                           Exposure, Encasement, &
                           MaterialData)
  Use dflib
  Use dfport
  Implicit None
!----------------------------------------------
! Arguments
!----------------------------------------------
  Character(30)  cProjName
  Character(30)  sec_name
  Character(120) cRunDir
  Integer(2)     Iter, Heat_mode, Protection
  Real(8)        rhoa, rHbeam, eps1res
  Real(8)        Timemax, Timestep, rpasimp
  Real(8)        Hbeam, Rbeam, AmsurV, Rdp
  Integer(2)     ProtRows
  Integer(2)     Analysis, Steel_grade, K_index, LoadType, Diagram
  Integer(2)     User, Rolled, Exposure, Encasement
  Real(8)        GammaMFi, NFid, Lfy, Lfz, K1_bend, K2_bend
  Real(8)        Lbraced, Psi, Zg, B_user, Tf_user
  Real(8)        Hi_user, Tw_user
  Real(8)        MaterialData(ProtRows,4)
!DEC$ ATTRIBUTES REFERENCE :: cProjName
!DEC$ ATTRIBUTES REFERENCE :: cRunDir
!DEC$ ATTRIBUTES REFERENCE :: sec_name
!DEC$ ATTRIBUTES REFERENCE :: MaterialData
!----------------------------------------------
! Local variables
!----------------------------------------------
  Integer       iochannel, iostat, ioresultL
  Character*43  Da2file, Da3file
  Character*30  cTemp
  Character*15  secnew_name
  Character*120 cTmpDir
  Character*8  Buffer
  Integer       i, j, i1, nbeleW
!
! CURRENT DATE
  Integer(2) tmpday, tmpmonth, tmpyear
!

  iochannel = 12

  cTmpDir = ""
  cTmpDir=cRunDir(1:index(cRunDir,char(00))-1)
  cTmpDir=Trim(cTmpDir)

  secnew_name = ""
  secnew_name = sec_name(1:index(sec_name,char(00))-1)
  secnew_name = Trim(secnew_name)

  cTemp = ""
  cTemp=cProjName(1:index(cProjName,char(00))-1)
  cTemp=Trim(cTemp)

! Number of steel iterations should be less than 30000.
! The Excel sheet can pass only 10 radial values at this point,
! but this is easy to change in the Excel VB macro. The idea is that
! for one room temperature one could calculate very many different
! steel calculations: different radial positions, different 
! protection, different sections, change some physical parameters
! of steel (density, steel grade,...).
! (If you change this, you must change it also in the CallOzone
!  subprogram)
  If ( Iter .le. 30000 ) Then
     Write (Buffer, '(i8)') Iter
     Da2file   = Trim(cTemp)//"-"//Trim(Adjustl(Buffer))
     Da3file   = Trim(Da2file)//".da3"
     Da2file   = Trim(Da2file)//".da2"
  Else
     Da2file   = Trim(cTemp)//"-ERROR"//".da2"
     Da3file   = Trim(cTemp)//"-ERROR"//".da3"
  End If
!
!----------------------------------------------
!  Change to the run directory
!----------------------------------------------
  ioresultL   = CHANGEDIRQQ (Trim(cTmpDir)) 
!
!----------------------------------------------
!  Create da2-file
!----------------------------------------------
!
  Open (iochannel, file = Trim(Da2file), action = 'write', form='formatted',iostat = iostat)
!
! CURRENT DATE: Check if this version is expired
!
  Call GETDAT(tmpyear, tmpmonth, tmpday)
  If ( (tmpyear.gt.2008) .or. ( (tmpmonth.gt.6) .and. (tmpyear.gt.2007) ) ) Then
     Write (iochannel,'(a,i2,''/'',i2.2,''/'',i4.4)') &
          'This version of PFS2 (CallOzone_Pfs_v2p1) has expired on ', &
          6, 30, 2008
     Stop
  End If

  Write (iochannel,'("GENERAL   ")')
  Write (iochannel,'("   Heating Protected")')
  Write (iochannel,'(2i10)') Heat_mode, Protection
  Write (iochannel,'("      rhoa        rH   epslres")')
  Write (iochannel,'(3f10.2)') rhoa, rHbeam, eps1res
  Write (iochannel,'("   Timemax  Timestep   rpasimp")')
  Write (iochannel,'(3f10.2)') Timemax, Timestep, rpasimp
  Write (iochannel,'("SECTION   ")')
  If ( User .ne. 1 ) Then
     Write (iochannel,'("      User")')
     Write (iochannel,'(i10)') User
  End If
  If ( Protection .eq. 0 ) Then
     If ( User .eq. 1 ) Then
        Write (iochannel,'("SectFactor")')
        Write (iochannel,'(1f10.2)') AmsurV
     Else
        Write (iochannel,'("      Name  Exposure")')
        Write (iochannel,'(a10,i10)') Trim(secnew_name), Exposure
     End If
  Else
!     Write (iochannel,'("    ApsurV       Rdp")')
!     Write (iochannel,'(2f10.2)') AmsurV, Rdp
     If ( User .eq. 1 ) Then
        Write (iochannel,'("SectFactor")')
        Write (iochannel,'(f10.2,f10.5)') AmsurV, Rdp
        Write (iochannel,'("      Rows")')
        Write (iochannel,'(i10)') ProtRows
     Else
        Write (iochannel,'("      Name  ExposureEncasement")')
        Write (iochannel,'(a10,2i10,f10.2)') Trim(secnew_name), Exposure, &
                                             Encasement
        Write (iochannel,'("PROTECTION")')
        Write (iochannel,'("       Rdp      Rows")')
        Write (iochannel,'(f10.5,i10)') Rdp, ProtRows
     End If
     Write (iochannel,'("    rTProV      rhop       rcp  rlambdap")')
     Do i = 1, ProtRows
        Write (iochannel,'(f10.2,f10.2,f10.2,f10.4)') (MaterialData(i,j),j = 1, 4)
     End Do
  End If
  Write (iochannel,'("LOCAL FIRE")')
  Write (iochannel,'("         H         r")')
  Write (iochannel,'(2f10.2)') Hbeam, Rbeam

  Close (iochannel)
!
!----------------------------------------------
!  Create da3-file
!----------------------------------------------
  Open (iochannel, file = Trim(Da3file), action = 'write', form='formatted',iostat = iostat)
  Write (iochannel,'("GENERAL   ")')
  Write (iochannel,'("  Analysis  GammaMFiSteelGrade")')
  Write (iochannel,'(i10,f10.2,i10)') Analysis, GammaMFi, Steel_grade
  Write (iochannel,'("LOAD      ")')
  If ( Analysis .ne. 3 ) Then
     Write (iochannel,'("      NFid")')
     Write (iochannel,'(f10.2)') NFid
  Else
     Write (iochannel,'("      MFid")')
     Write (iochannel,'(f10.2)') NFid
  End If
  Write (iochannel,'("SECTION   ")')
  Write (iochannel,'("      User")')
  Write (iochannel,'(i10)') User
  If ( User .eq. 1 ) Then
     Write (iochannel,'("         B        Tf        Hi        Tw        Tw")')
     Write (iochannel,'(4f10.2,i10)') B_user, Tf_user, Hi_user, Tw_user, Rolled
  Else
     Write (iochannel,'("      Name")')
     Write (iochannel,'(a10)') Trim(secnew_name)
  End If

  If ( Analysis .eq. 2 ) Then
     Write (iochannel,'("COMPRESSION")')
     Write (iochannel,'("       Lfy       Lfz")')
     Write (iochannel,'(2f10.2)') Lfy, Lfz
  End If

  If ( Analysis .eq. 3 ) Then
     Write (iochannel,'("BENDING   ")')
     Write (iochannel,'("   K index        K1        K2   Lbraced")')
     Write (iochannel,'(i10,3f10.2)') K_index, K1_bend, K2_bend, Lbraced
     Write (iochannel,'("  LoadType")')
     Write (iochannel,'(i10)') LoadType
     If ( LoadType .eq. 1 ) Then
        Write (iochannel,'("       Psi")')
        Write (iochannel,'(f10.2)') Psi
     Else
        Write (iochannel,'("   Diagram        Zg")')
        Write (iochannel,'(i10,f10.2)') Diagram, Zg
     End If
  End If

  Close (iochannel)
!--------------------------------------------------
! Finished so return
!--------------------------------------------------
  Return
End Subroutine WriteSteelData
!***********************************************************
!
!  Read Ozone results: wall temperatures, openings, transition times
!
!  2nd December 2002
!  This works if:
!     -  exactly four walls
!     -  each partition has 7 elements
!     -  max. number of vents is 10
!     -  max. number of materials is 10
!
!DEC$ ATTRIBUTES STDCALL,DLLEXPORT        :: ReadWallEtc
!DEC$ ATTRIBUTES ALIAS : "ReadWallEtc"   :: ReadWallEtc
!***********************************************************
Subroutine ReadWallEtc (cProjName, cRunDir, cProjectBase, CallModeIn, &
     iStrategy,  &
     TFlashOver, rTFuelIgni, rPercAfionAf, rPercZsonH, &
     rAfloor, rHeight, zFireElev, rfuelheight, &
     Tini, rPinitial, rDHc, repsilglass, timemax, &
     iceiling, ifloor, nbwall, nbmat, nbopen, irmf, &
     nbptrmf, iceiltype, icomb, iairent, iExtPrint, &
     ropper1, ropT2, ropper2, ropT3, ropper3, &
     roptime2, ropT0, openings, materials, MaxWallNodes, &
     nbeleC, nbeleF, walldata, partdata, rfirearea, QcriMJ, Tinscri, rCrit, &
     OpData )

  Use dflib
  Use dfport
  Implicit None
!----------------------------------------------
! Arguments
!----------------------------------------------
  Character(30)  cProjName, cProjectBase
  Character(120) cRunDir
  Integer(2)     iStrategy, CallModeIn
  Real(8)        TFlashOver, rTFuelIgni, rPercAfionAf, rPercZsonH  
  Real(8)        rAfloor, rHeight, zFireElev, rfuelheight
  Real(8)        Tini, rPinitial, rDHc, repsilglass
  Real(8)        Timemax
  Integer(2)     iceiling, ifloor, nbwall, nbmat, nbopen, irmf
  Integer(2)     nbptrmf, iceiltype, icomb, iairent, iExtPrint
  Real(8)        ropper1(3), ropT2(2), ropper2(2), ropT3(2), ropper3(2)
  Real(8)        roptime2, ropT0
  Real(8)        openings(10,7), materials(10,7)
  Integer(2)     MaxWallNodes, nbeleC, nbeleF
  Real(8)        rfirearea
  Real(8)        walldata(2*nbwall)
!  Real(8)        partdata(7,2*nbwall+4)   ! Old one, now more than 7 lines
  Real(8)        partdata(MaxWallNodes,2*nbwall+4)
  Real(8)        Tinscri(6)
  Real(8)        rCrit(4), OpData(2,nbopen), QcriMJ
!DEC$ ATTRIBUTES REFERENCE :: cProjName, cRunDir, cProjectBase
!DEC$ ATTRIBUTES REFERENCE :: ropper1, ropT2, ropper2, ropT3, ropper3
!DEC$ ATTRIBUTES REFERENCE :: openings, materials
!DEC$ ATTRIBUTES REFERENCE :: walldata, partdata, tinscri, rCrit, OpData
!----------------------------------------------
! Local variables
!----------------------------------------------
  Integer       iochannel, iostat, ioresultL, ioresultI, dbgchannel, islabchannel
  Character*34  DbgFile, CriFile, Wu1File, Wu2File, Wu3File, Wu4File, CeiFile, OpFile
  Character*50  cTmpSave
  Character*80  line
  Character*30  cTemp
  Character*120 cTmpDir
  Character*14  criteria
  Logical Debug, Slab
  Integer       i, j, i1, nbeleW, iLines, iwall, nele
  Integer       iop, ivar, CallMode, iSlabCalc
  Real*8        tempcri, timecri, xcrit, tcrit
  Real*8        epsi, rcon, tinsu
  Real*8, Allocatable :: InData(:,:), tgas(:), t(:)
!
! CURRENT DATE
  Integer(2) tmpday, tmpmonth, tmpyear
!
  CallMode = CallModeIn
  iochannel    = 22
  islabchannel = 23
  dbgchannel   = 25

  cTmpDir = ""
  cTmpDir=cRunDir(1:index(cRunDir,char(00))-1)
  cTmpDir=Trim(cTmpDir)

  cTmpSave = ""
  cTmpSave=cProjectBase(1:index(cProjectBase,char(00))-1)
  cTmpSave=Trim(cTmpSave)

  cTemp = ""
  cTemp=cProjName(1:index(cProjName,char(00))-1)
  cTemp=Trim(cTemp)
  DbgFile   = Trim(cTmpSave)//".dbg"
!!!  DbgFile   = Trim(cTemp)//".dbg"
  CriFile   = Trim(cTemp)//".cri"
  OpFile    = Trim(cTemp)//".op"
  CeiFile   = Trim(cTemp)//".cei"
  Wu1File   = Trim(cTemp)//".WU1"
  Wu2File   = Trim(cTemp)//".WU2"
  Wu3File   = Trim(cTemp)//".WU3"
  Wu4File   = Trim(cTemp)//".WU4"
  rCrit = 0.d0
  Tinscri = 0.d0
  OpData = 0.d0

!
!----------------------------------------------
!  Change to the run directory
!----------------------------------------------
  ioresultL   = CHANGEDIRQQ (Trim(cTmpDir)) 


  iSlabCalc = CallMode/10
  iSlabCalc = Min(1,iSlabCalc)    ! Should be 0 or 1
! iSlabCalc 0: Normal steel beam resistance calculation
!           1: BRE composite slab method (Bailey)
  CallMode  = Mod(CallMode,10)    ! Should be 0 - 5
  If ( CallMode .lt. 0 ) CallMode = 0
  If ( CallMode .gt. 6 ) CallMode = 6
! CallMode 0: gas temperature
!          1: gas temperature, debug mode
!          2: gas and steel (and slab) temperatures
!          3: gas and steel (and slab) temperatures, debug mode
!          4: gas temperatures, save, no debug
!          5: gas and steel (and slab) temperatures, save, no debug
!          6: as 5 but debug mode
  If ( CallMode .ge. 4 ) Then
     If ( CallMode .eq. 4 ) CallMode = 0
     If ( CallMode .eq. 5 ) CallMode = 2
     If ( CallMode .ge. 6 ) CallMode = 3
  End If

  ! No composite slab calculation, if only gas temperature is calculated
  If ( (CallMode.lt.2) .or. (Callmode.eq.4) ) iSlabCalc = 0

! CallMode 0: gas temperature
!          1: gas temperature, debug mode
!          2: gas and steel temperatures
!          3: gas and steel temperatures, debug mode
  If ((CallMode .eq. 1) .or. (CallMode .eq. 3)) Then
     Debug = .true.
  Else
     Debug = .false.
     ioresultI = DELFILESQQ (Trim(cTemp)//".WL1") 
     ioresultI = DELFILESQQ (Trim(cTemp)//".WL2") 
     ioresultI = DELFILESQQ (Trim(cTemp)//".WL3") 
     ioresultI = DELFILESQQ (Trim(cTemp)//".WL4") 
     ioresultI = DELFILESQQ (Trim(cTemp)//".flo") 
  End If

  If ( iSlabCalc .eq. 1 ) Then
     Slab = .true.
  Else
     Slab = .false.
  End If

  If (Debug) Then
     Open (dbgchannel, file=Trim(DbgFile), action='readwrite', position='append', form='formatted', iostat=iostat)
     Write (dbgchannel,*) '******* ReadWallEtc *******'
     Write (dbgchannel,*) 'Change dir, ioresult is ', ioresultL
     Write (dbgchannel,*) 'cTmpDir: ', cTmpDir
  End If
!
! CURRENT DATE: Check if this version is expired
!
  Call GETDAT(tmpyear, tmpmonth, tmpday)
  If ( (tmpyear.gt.2008) .or. ( (tmpmonth.gt.6) .and. (tmpyear.gt.2007) ) ) Then
     If (Debug) Then
        Write (dbgchannel,*)
        Write (dbgchannel,'(a,i2,''/'',i2.2,''/'',i4.4)') &
             'This version of PFS2 (CallOzone_Pfs_v2p1) has expired on ', &
             6, 30, 2008
        Write (dbgchannel,*)
     End If
     Stop
  End If

!
!----------------------------------------------
!  Open transition criteria file
!----------------------------------------------
  Open (iochannel, file = Trim(CriFile), action = 'read', form='formatted',iostat = iostat)
  Read (iochannel, '(A)') line
  Do While (.not. EOF(iochannel))
     Read (iochannel, fmt='(a,2f10.2)') criteria, tempcri, timecri
     Select Case (Trim(ADJUSTL(criteria)))
     Case ('CRITERION C1:')
        rCrit (1) = timecri   ! flashover temperature reached
     Case ('CRITERION C2:')
        rCrit (2) = timecri   ! all fuel is burning ( T > Tign )
     Case ('CRITERION C3:')  ! 2zone to 1zone, hot layer too close to floor
        rCrit (3) = timecri
     Case ('CRITERION C4:')  ! 2zone to 1zone, fuel area too large
        rCrit (4) = timecri
     Case Default
        Continue
     End Select
     If (Debug) Then
        Write (dbgchannel,*) 'criteria, tempcri, timecri'
        Write (dbgchannel,*) criteria, tempcri, timecri
     End If
  End Do
  Close (iochannel)
  If (Debug) Then
     Write (dbgchannel,*) ' cri file closed'
  Else
     ioresultI = DELFILESQQ (CriFile) 
  End If

!----------------------------------------------
! Wall and ceiling: insulation criteria
!----------------------------------------------
! Floor file is not opened, only walls and ceiling are considered

!----------------------------------------------
! Loop over the wall data files 
!----------------------------------------------
  iwall = 1
  If ( nbwall .ne. 4 ) Then
!     Write (25,*) 'Error in ReadWallEtc: nbwall ne 4 ', nbwall
     Stop
  End If
  Do iwall = 1, nbwall
     Select Case (iwall)
     Case (1)
        Open (iochannel, file = Trim(Wu1File), action = 'read', form='formatted',iostat = iostat)
     Case (2)
        Open (iochannel, file = Trim(Wu2File), action = 'read', form='formatted',iostat = iostat)
     Case (3)
        Open (iochannel, file = Trim(Wu3File), action = 'read', form='formatted',iostat = iostat)
     Case (4)
        Open (iochannel, file = Trim(Wu4File), action = 'read', form='formatted',iostat = iostat)
     Case Default
        If (Debug) Then
           Write (dbgchannel,*) 'Error, ReadWallEtc, wall file'
        End If
        Stop
     End Select
     i = 0
     Do While (.not. EOF(iochannel))
        i = i + 1
        Read (iochannel, *)
     End Do
     ! First line is the title line
     iLines = i - 1
     Allocate( InData(iLines,Nint(walldata((iwall-1)*2+2))+1) )
     Allocate( t(iLines) )
     Allocate( tgas(iLines) )
     Rewind (iochannel)
     
     Read (iochannel, '(A)') line
     Do i = 1, iLines
        Read (iochannel, fmt='(f8.1,9f8.2)') t(i), tgas(i),   &
             (InData(i,j), j = 1, Nint(walldata((iwall-1)*2+2))+1 )
     End Do
     nele = Nint(walldata((iwall-1)*2+2))
     epsi = materials(Nint(partdata(1,(iwall-1)*2+6)),4)
     rcon = materials(Nint(partdata(1,(iwall-1)*2+6)),6)
     Call Insulation(Timemax, Tini, epsi, rcon, iLines, nele, QcriMJ, tgas, &
          t, InData, tinsu )
! First is floor, then ceiling, then the walls
     Tinscri(iwall+2) = tinsu
     If (Debug) Then
        Write (dbgchannel,*) 'Wall: insulation failure (s): ', tinsu
     End If
     Deallocate( tgas )
     Deallocate( t )
     Deallocate( InData )
     Close (iochannel)
     If (Debug) Then
        Write (dbgchannel,*) ' wall file closed, wall no ', iwall
     Else
        Select Case (iwall)
        Case (1)
           ioresultI = DELFILESQQ (Wu1File) 
        Case (2)
           ioresultI = DELFILESQQ (Wu2File) 
        Case (3)
           ioresultI = DELFILESQQ (Wu3File) 
        Case (4)
           ioresultI = DELFILESQQ (Wu4File) 
        Case Default
           If (Debug) Then
              Write (dbgchannel,*) 'Error2, ReadWallEtc, wall file'
           End If
        End Select
     End If
  End Do   ! loop over the (upper) walls

!----------------------------------------------
! Open the ceiling data file (Floor is not treated here)
!----------------------------------------------
  Open (iochannel, file = Trim(CeiFile), action = 'read', form='formatted',iostat = iostat)
  i = 0
  Do While (.not. EOF(iochannel))
     i = i + 1
     Read (iochannel, *)
  End Do
  ! First line is the title line
  iLines = i - 1
  Allocate( InData(iLines,nbeleC+1) )
  Allocate( t(iLines) )
  Allocate( tgas(iLines) )
  Rewind (iochannel)
     
  Read (iochannel, '(A)') line
  tgas=0.d0
  ! If debug-mode and slab-mode and MC simulation then one should make a 
  ! copy of the  _imc.cei file (if no debug, rename is enough). A copy is
  ! made by write statements, because no copyfileqq (or similar) fortran
  ! commands exists.
  If ( Debug .and. Slab ) Then
     If ( .not.( Trim(cTmpSave) .eq. Trim(cTemp) ) ) Then
        ioresultI = DELFILESQQ (Trim(cTmpSave)//".cei") 
        Open (islabchannel, file = Trim(cTmpSave)//".cei", &
             action = 'write', form='formatted',iostat = iostat)
        Write (islabchannel, '(A)') line
     End If
  End If

  Do i = 1, iLines
     Read (iochannel, fmt='(f8.1,9f8.2)') t(i), tgas(i),   &
          (InData(i,j), j = 1, nbeleC+1 )

     If ( Debug .and. Slab ) Then
        If ( .not.( Trim(cTmpSave) .eq. Trim(cTemp) ) ) Then
           Write (islabchannel, fmt='(f8.1,9f8.2)') t(i), tgas(i),   &
                (InData(i,j), j = 1, nbeleC+1 )
        End If
     End If

  End Do

  If ( Debug .and. Slab ) Then
     If ( .not.( Trim(cTmpSave) .eq. Trim(cTemp) ) ) Then
        Close (islabchannel)
     End If
  End If

  nele = nbeleC
  epsi = materials(Nint(partdata(1,2)),4)
  rcon = materials(Nint(partdata(1,2)),6)
  Call Insulation(Timemax, Tini, epsi, rcon, iLines, nele, QcriMJ, tgas, &
       t, InData, tinsu )
  ! First is ceiling, then floor, then the walls
  Tinscri(1) = tinsu
  Tinscri(2) = 0.d0   ! floor is not calculated
  If (Debug) Then
     Write (dbgchannel,*) 'Ceiling: insulation failure (s): ', tinsu
  End If
!
  If ( rfirearea .ge. rPercAfionAf*rAfloor ) Then
! ozonepfs.dll BUG FIX: non-localized and T > Tflashover  (C1)
! This happens, when T> Tfo is reached first and then the calculation
! switches to one zone and all fuel is burning, so there can not be
! any additional transitions anymore  ==> All rCrit(1-4) are zero
! when this bug is happening.
     If (Debug) Then
        Write (dbgchannel,*) 'rfirearea    ', rfirearea
        Write (dbgchannel,*) 'rAfloor      ', rAfloor
        Write (dbgchannel,*) 'rPercAfionAf ',rPercAfionAf
     End If
     If ( ( rCrit(1).le.1.d-2 .and. rCrit(2).le.1.d-2 .and.  &
          rCrit(3).le.1.d-2 .and. rCrit(4).le.1.d-2      ) &
          .and. (Maxval(tgas) .ge. TFlashOver) ) Then
        ! Find the flashover time
        Do i = 1, iLines - 1
           If ( tgas(i+1) .ge. TFlashOver ) Then
              If ( tgas(i+1)-tgas(i) .gt. 1.d-3 ) Then
                 rCrit(1) = t(i) + (t(i+1)-t(i)) * &
                      (TFlashOver-tgas(i)) / &
                      (tgas(i+1)-tgas(i))
              Else
                 rCrit(1) = t(i)   ! do not divide with zero
              End If
              If ( rCrit(4) .le. 1.d-2 ) rCrit(4) = rCrit(1)
              Exit
           End If
        End Do
     End If
     ! If all fuel starts to burn then also fire area is larger then
     ! the trasition 2zone to 1zone criterion (distributed loads)
     If ( rCrit(2).gt.1.d-2 .and. rCrit(4).le.1.d-2) rCrit(4) = rCrit(2)

  End If  ! distributed fire load: bug fix

  Deallocate( tgas )
  Deallocate( t )
  Deallocate( InData )
  Close (iochannel)

  If (Debug) Then
     Write (dbgchannel,*) ' Ceiling file closed'

  Else

     If ( Slab ) Then
        If ( .not.( Trim(cTmpSave) .eq. Trim(cTemp) ) ) Then
           ! a MC iteration, i.e. casename = projectname_imc
           ! Interactive run: casename = projectname
           ioresultI = DELFILESQQ (Trim(cTmpSave)//".cei") 
           ioresultI = RENAMEFILEQQ(Trim(CeiFile), Trim(cTmpSave)//".cei")
        End If
     Else
        ioresultI = DELFILESQQ (Trim(CeiFile)) 
     End If

  End If    ! Debug
!----------------------------------------------
! Wall and ceiling files are done
!----------------------------------------------

!----------------------------------------------
! Window opening times
!----------------------------------------------
  OpData = 9.d9
  Open (iochannel, file = Trim(OpFile), action = 'read', form='formatted',iostat = iostat)
  Read (iochannel, '(A)') line
  If (Debug) Then
     Write (dbgchannel,*) 'Opening file opened'
  End If
  Do While (.not. EOF(iochannel))
     Read (iochannel, fmt='(2i5,2f10.2)') iop, ivar, xcrit, tcrit
     If (Debug) Then
        Write (dbgchannel, fmt='(2i5,2f10.2)') iop, ivar, xcrit, tcrit
     End If
! There is at most two criteria for each window. Note, that Ozone
! output may have many opening times for each criteria, e.g. if the
! rises first and the drops below the temperature criteria. (Program
! writes times every time that the temperature value is crossed.)
! 'Bug fix': Only the first time is needed, because after that the
! window is already broken. Note: Criterion 3 (time dependent) is not
! written to the file .op, because the time is predetermined in the
! input.
     If ( OpData(1,iop) .gt. 1.d9 ) Then
        OpData(1,iop) = tcrit
     Else If ( OpData(2,iop) .gt. 1.d9 ) Then
        OpData(2,iop) = tcrit
     End If
     If (Debug) Then
        Write (dbgchannel,*) iop,OpData(1,iop),OpData(2,iop)
     End If
  End Do
  Do iop = 1, nbopen
     If ( openings(iop,6) .eq. 3 ) Then
        OpData(1,iop) = roptime2
        OpData(2,iop) = 0.d0
     End If
     If ( OpData(1,iop) .ge. 1.d9 ) Then
        OpData(1,iop) = 0.d0
     End If
     If ( OpData(2,iop) .ge. 1.d9 ) Then
        OpData(2,iop) = 0.d0
     End If
  End Do
  Close (iochannel)
  If (Debug) Then
     Write (dbgchannel,*) ' Opening file closed'
  Else
     ioresultI = DELFILESQQ (OpFile) 
  End If

! Format of the .op file:
! Opening number, type of variation, criterion value, time of happening
!
!--------------------------------------------------
! Finished so clean up and return
!--------------------------------------------------
  If (Debug) Then
     Write (dbgchannel,*) ' Debug file closed'
     Write (dbgchannel,*) '****** ReadWallEtc ends ********'
  End If
  Close (dbgchannel)
  Return

Contains

  Subroutine Insulation(Timemax, Tini, epsi, rcon, iLines, nele, QcriMJ, tgas, &
                        t, InData, tinsu)
! This subroutine calculates the time when the critical heat flux is 
! reached (QcriMJ). Both convective and radiative heating are treated.

! Passed varialbles
    Integer iLines, nele
    Real*8 Timemax, Tini, epsi, rcon, QcriMJ
    Real*8 tgas(iLines), t(iLines), InData(iLines,nele)

! Local variables
    Real*8 sigma, sum, qdotin, Tabs, tinsu, criticalT

    sigma = 5.6704d-8   ! Stefan-Bolzmann constant
    Tabs  = 273.15d0    ! Celsius to Kelvin

    If ( QcriMJ .gt. 0 ) Then
! Trapezoidal rule for the total heat flux into the wall
       sum = 0.d0
       tinsu = Timemax
!         tinsu = t(iLines)
       Do i = 1, iLines-1
          qdotin = epsi*sigma*((tgas(i+1)+Tabs)**4 - (InData(i+1,1)+Tabs)**4) + &
               rcon*(tgas(i+1) - InData(i+1,1)) + &
               epsi*sigma*((tgas(i)+Tabs)**4 - (InData(i,1)+Tabs)**4) + &
               rcon*(tgas(i) - InData(i,1))
          qdotin = 0.5d0*qdotin  ! Qdot at the mid-point
! If the criteria is reached during this time interval, interpolate
! the failure time and exit the do-loop.
          If ( sum + qdotin*(t(i+1)-t(i)) .ge. QcriMJ*1.0d6 ) Then
             tinsu = t(i) +  ( QcriMJ*1.0d6 - sum ) / ( qdotin ) 
             Exit
          End If
          sum = sum + qdotin*( t(i+1)-t(i) )
       End Do
    Else
! Temperature criterion for the back side node (Tcrit = -QcriMJ)
       criticalT = - QcriMJ
       tinsu = Timemax
!         tinsu = t(iLines)
       Do i = 1, iLines - 1
          If ( InData(i+1,nele) .ge. criticalT ) Then
             If ( InData(i+1,nele)-InData(i,nele) .gt. 1.d-3 ) Then
                tinsu = t(i) + (t(i+1)-t(i)) * &
                     (criticalT-InData(i,nele)) / &
                     (InData(i+1,nele)-InData(i,nele))
             Else
                tinsu = t(i)   ! do not divide with zero
             End If
             Exit
          End If
       End Do
    End If
    Return
  End Subroutine Insulation

End Subroutine ReadWallEtc
!***********************************************************

!***********************************************************
!
!  Read OZone Dat -file
!
!  4th November 2003
!  This works if:
!     -  t2-fire is used (user given fire is not read in, all the
!        other parameters are read in)
!     -  exactly four walls
!     -  each partition has <= PMaxWallNodes elements
!     -  max. number of vents is 10
!     -  max. number of materials is 10
!
!DEC$ ATTRIBUTES STDCALL,DLLEXPORT        :: ReadOzoneData
!DEC$ ATTRIBUTES ALIAS : "ReadOzoneData"   :: ReadOzoneData
!***********************************************************
Subroutine ReadOzoneData (cProjName, cRunDir, iStrategy,  &
     TFlashOver, rTFuelIgni, rPercAfionAf, rPercZsonH, &
     rAfloor, rHeight, cHeight, zFireElev, rfuelheight, &
     Tini, rPinitial, rDHc, repsilglass, &
     xtimestep, printstep, timemax, &
     iceiling, ifloor, nbwall, nbmat, nbopen, irmf, &
     nbptrmf, iceiltype, icomb, iairent, nbsmokeext, &
     nbHVgroup, iExtPrint, &
     ropper1, ropT2, ropper2, ropT3, ropper3, &
     roptime2, ropT0, openings, materials, fvents, hvents, &
     MaxWallNodes, nbeleC, nbeleF, walldata, partdata, rmfmax, tplateau, &
     tendplate, tendfire, rfirearea, RHRf, Pnbwall, PMaxWallNodes )

  Use dflib
  Use dfport
  Implicit None
  Integer(2) Pnbwall, PMaxWallNodes
!  Parameter ( Pnbwall = 4 )   ! Exactly four walls
!  Parameter ( PMaxWallNodes = 50 )   ! Maximum input wall nodes
!----------------------------------------------
! Arguments
!----------------------------------------------
  Character(30)  cProjName
  Character(120) cRunDir
  Integer(2)     iStrategy
  Real(8)        TFlashOver, rTFuelIgni, rPercAfionAf, rPercZsonH  
  Real(8)        rAfloor, rHeight, cHeight, zFireElev, rfuelheight
  Real(8)        Tini, rPinitial, rDHc, repsilglass
  Real(8)        xtimestep, printstep, timemax
  Integer(2)     iceiling, ifloor, nbwall, nbmat, nbopen, irmf
  Integer(2)     nbptrmf, iceiltype, icomb, iairent, nbsmokeext
  Integer(2)     nbHVgroup, iExtPrint
  Real(8)        ropper1(3), ropT2(2), ropper2(2), ropT3(2), ropper3(2)
  Real(8)        roptime2, ropT0
  Real(8)        openings(10,7), materials(10,7), fvents(3,5), hvents(3,4)
  Integer(2)     MaxWallNodes, nbeleC, nbeleF
  Real(8)        walldata(2*Pnbwall)
!  Real(8)        partdata(7,2*nbwall+4)   ! This is old, now more than 7 lines
  Real(8)        partdata(PMaxWallNodes,2*Pnbwall+4)
  Real(8)        rmfmax, tplateau, tendplate, tendfire, rfirearea, RHRf
!DEC$ ATTRIBUTES REFERENCE :: cProjName, cRunDir, iStrategy
!DEC$ ATTRIBUTES REFERENCE :: TFlashOver, rTFuelIgni, rPercAfionAf, rPercZsonH
!DEC$ ATTRIBUTES REFERENCE :: rAfloor, rHeight, cHeight, zFireElev, rfuelheight
!DEC$ ATTRIBUTES REFERENCE :: Tini, rPinitial, rDHc, repsilglass
!DEC$ ATTRIBUTES REFERENCE :: xtimestep, printstep, timemax
!DEC$ ATTRIBUTES REFERENCE :: iceiling, ifloor, nbwall, nbmat, nbopen, irmf
!DEC$ ATTRIBUTES REFERENCE :: nbptrmf, iceiltype, icomb, iairent, nbsmokeext
!DEC$ ATTRIBUTES REFERENCE :: nbHVgroup
!DEC$ ATTRIBUTES REFERENCE :: ropper1, ropT2, ropper2, ropT3, ropper3
!DEC$ ATTRIBUTES REFERENCE :: roptime2, ropT0, openings, materials, fvents, hvents
!DEC$ ATTRIBUTES REFERENCE :: MaxWallNodes, nbeleC, nbeleF
!DEC$ ATTRIBUTES REFERENCE :: walldata, partdata, rmfmax, tplateau
!DEC$ ATTRIBUTES REFERENCE :: tendplate, tendfire, rfirearea
!DEC$ ATTRIBUTES REFERENCE :: RHRf
!----------------------------------------------
! Local variables
!----------------------------------------------
  Integer       iochannel, iostat, ioresultL
  Character*34  DatFile
  Character*30  cTemp
  Character*120 cTmpDir
  Integer       i, j, i1, nbeleW
  Character*160 Buffer
  !
  ! CURRENT DATE
  Integer(2) tmpday, tmpmonth, tmpyear
  !

  openings     = 0.d0
  materials    = 0.d0
  fvents       = 0.d0
  hvents       = 0.d0
  walldata     = 0.d0
  partdata     = 0.d0
  MaxWallNodes = 0
  iochannel = 12
  cTmpDir = ""
  cTmpDir=cRunDir(1:index(cRunDir,char(00))-1)
  cTmpDir=Trim(cTmpDir)
  cTemp = ""
  cTemp=cProjName(1:index(cProjName,char(00))-1)
  cTemp=Trim(cTemp)
  DatFile   = Trim(cTemp)
!  DatFile   = Trim(cTemp)//".dat"
  !
  !----------------------------------------------
  !  Change to the run directory
  !----------------------------------------------
  ioresultL   = CHANGEDIRQQ (Trim(cTmpDir)) 
  !
  !----------------------------------------------
  !  Open datfile
  !----------------------------------------------
  Open (iochannel, file = Trim(DatFile), action = 'read', form='formatted',iostat = iostat)
  !
  ! CURRENT DATE: Check if this version is expired
  !
  Call GETDAT(tmpyear, tmpmonth, tmpday)

  ! The name of the project is not read in
  Read (iochannel,*)
  !  cTemp = Trim(cTemp)
  !  Write (6,'(" iStrategy")')
  Read (iochannel,*)
  Read (iochannel,*) iStrategy

  !  Write (6,'("     TFlashOver     rTFuelIgni   rPercAfionAf     rPercZsonH ")')
  Read (iochannel,*)
  Read (iochannel,*) TFlashOver, rTFuelIgni, rPercAfionAf, rPercZsonH

  !  Write (6,'("   rAfloor        rHeight      zFireElev    rfuelheight ")')
  Read (iochannel,*)
  Read (iochannel,*) rAfloor, rHeight, zFireElev, rfuelheight

  !  Write (6,'("      Tini      rPinitial           rDHc    repsilglass ")')
  Read (iochannel,*)
  Read (iochannel,*) Tini, rPinitial, rDHc, repsilglass

  !  Write (6,'(" xtimestep      printstep        timemax ")')
  Read (iochannel,*)
  Read (iochannel,*) xtimestep, printstep, timemax

  ! Switches and numbers of different objects
  !  Write (6,'(a)') '  iceiling    ifloor    nbwall     nbmat    nbopen' // &
  !       '      irmf   nbptrmf iceiltype     icomb   iairent ' // &
  !       'nbsmokeext nbHVgroup iExtPrint '
  Read (iochannel,*)
  Read (iochannel,'(13i10)') iceiling, ifloor, nbwall, nbmat, nbopen, irmf, nbptrmf, &
          iceiltype, icomb, iairent, nbsmokeext, nbHVgroup, iExtPrint
  If ( nbmat .gt. 10 ) Then
     Stop 'nbmat is greater than 10'
  End If
  If ( nbopen .gt. 10 ) Then
     Stop 'nbopen is greater than 10'
  End If
  If ( nbwall .ne. 4 ) Then
     Stop 'nbwall is not equal to 4'
  End If

  ! Finite element data for the partitions
  !  Write (6,'(a)') '   CEILING'
  !  Write (6,'(a)') '    nbeleC'
  !  Write (6,'(a)') '    rlelem   nmat'
  Read (iochannel,*)
  Read (iochannel,*)
  Read (iochannel,*)

  If (iceiltype .eq. 0) Then
     Read (iochannel,*) nbeleC
  Else
     Read (iochannel,*) nbeleC, cHeight
  End If

  If ( nbeleC .gt. PMaxWallNodes ) Then
     Stop 'Too many ceiling nodes, PMaxWallNodes'
  End If
  MaxWallNodes = Max( MaxWallNodes, nbeleC )

  Do j = 1, nbeleC
     Read (iochannel,*) partdata(j,1), partdata(j,2)
  End Do

  !  Write (6,'(a)') '     FLOOR'
  !  Write (6,'(a)') '    nbeleF'
  !  Write (6,'(a)') '    rlelem   nmat'
  Read (iochannel,*)
  Read (iochannel,*)
  Read (iochannel,*)

  If (iceiltype .eq. 0 .or. iceiltype .eq. 2) Then
     Read (iochannel,*) nbeleF
  Else
     Read (iochannel,*) nbeleF, cHeight
  End If

  If ( nbeleF .gt. PMaxWallNodes ) Then
     Stop 'Too many floor nodes, PMaxWallNodes'
  End If
  MaxWallNodes = Max( MaxWallNodes, nbeleF )

  Do j = 1, nbeleF
     Read (iochannel,*) partdata(j,3), partdata(j,4)
  End Do

  !  Write (6,'(a)') '     WALLS'
  !  Write (6,'(a)') 'walllength nbeleW'
  !  Write (6,'(a)') '    rlelem   nmat'
  Read (iochannel,*)
  Read (iochannel,*)
  Read (iochannel,*)

  Do i = 1, nbwall
     i1 = 2*(i-1) + 1 
     Read (iochannel,*) walldata(i1), walldata(i1+1)
     nbeleW = Nint(walldata(i1+1))
     If ( nbeleW .gt. PMaxWallNodes ) Then
        Stop 'Too many wall nodes, PMaxWallNodes'
     End If
     MaxWallNodes = Max( MaxWallNodes, nbeleW )

     i1 = i1 + 4 ! now the first 4 points are used for floor and ceiling
     Do j = 1, nbeleW
        Read (iochannel,*) partdata(j,i1), partdata(j,i1+1)
     End Do
  End Do

  ! Materials used for the walls, floor and ceiling
  !  Write (6,'(a)') 'MATERIAL FOR PARTITIONS'
  Read (iochannel,*)

  !  Write (6,'(a)') '       rho        rk        rc repsilint' // &
  !       ' repsilext     rhint     rhext'
  Read (iochannel,*)

  nbmat = 0
  Do
     Read (iochannel, '(A)') Buffer
     If ( Buffer(1:10) .eq. '     rsill' ) Then
        Exit
     End If
     nbmat = nbmat + 1
     Read (Buffer,*) (materials(nbmat,j), j=1,7)
  End Do

  ! Openings (doors and windows)
  !  Write (6,'("     rsill   rsoffit    rwidth    rcvent   nopwall   iopevol     iadia ")')
  !  Read (iochannel,*)

  nbopen = 0
  Do
     Read (iochannel, '(A)') Buffer
     If ( (Buffer(1:10) .eq. '    rmfmax') .or. (Buffer(1:10) .eq. '      time') ) Then
        Exit
     End If
     nbopen = nbopen + 1
     Read (Buffer,*) (openings(nbopen,j), j=1,7)
  End Do

  ! RHR data
  If ( irmf .eq. 0 ) Then
     !     Write (6,'("    rmfmax  tplateau tendplate  tendfire rfirearea      RHRf ")')
     !     Read (iochannel,*)

     Read (iochannel,*) rmfmax, tplateau, tendplate, &
          tendfire, rfirearea, RHRf
     RHRf = RHRf/1000.d0
     Read (iochannel,*)
  Else
     !
     ! Now just skip the user given time series.
     Do
        Read (iochannel, '(A)') Buffer
        If ( Buffer(1:10) .eq. '   ropper1' ) Then
           Exit
        End If
     End Do

  End If

  ! Openings: Opening criteria:
  ! first line:  temperature dependent, stairs
  ! second line: temperature dependent, piecewise linear
  ! third line:  time dependent
  ! fourth line: temperature dependent, single step, i.e. closed vs open

  !  Write (6,'("   ropper1     ropT2   ropper2     ropT3   ropper3 ")')
  Read (iochannel,*) ropper1(1), ropT2(1), ropper2(1), ropT3(1), ropper3(1)

  Read (iochannel,*) ropper1(2), ropT2(2), ropper2(2), ropT3(2), ropper3(2)
  !  Write (6,'("   ropper1  roptime2 ")')
  Read (iochannel,*)
  Read (iochannel,*) ropper1(3), roptime2

  !  Write (6,'("               ropT0 ")')
  Read (iochannel,*)
  Read (iochannel,*) ropT0

  ! Forced vents (vertical)
!  Write (6,'("       zFV    DiamFV       VFV   signMFV iopevolFV ")')
  Read (iochannel,*)
  Do i = 1, nbsmokeext
     Read (iochannel,*) (fvents(i,j), j=1,5)
  End Do

  ! Horizontal vents (smoke hatches)
  !  Write (6,'("      nbHV   rAreaHV iopevolHV ")')
  Read (iochannel,*)
  Do i = 1, nbHVgroup
     Read (iochannel,*) (hvents(i,j), j=1,4)
  End Do
!  Read (iochannel,*)
!--------------------------------------------------
! Ready to finish
!--------------------------------------------------
  Close (iochannel)
!--------------------------------------------------
! Finished so return
!--------------------------------------------------
  Return
End Subroutine ReadOzoneData
!***********************************************************
